COMMENT ⊗   VALID 00005 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	begin "unplot" comment An algorithm to undiscretize data.
C00006 00003	Simplex procedure
C00013 00004	The main algorithm
C00017 00005	The driver program
C00019 ENDMK
C⊗;
begin "unplot" comment An algorithm to undiscretize data.

Given points $a↓0$, $\ldotss$, $a↓{n+1}$ (typically all integers with
$|a↓{j+1}-a↓j|≤1$, but this is not required), this algorithm finds
points $z↓0$, $\ldotss$, $z↓{n+1}$ such that $|z↓j-a↓j|≤{1\over2}$ for all
$j$, and also $z↓j=a↓j$ for $j = 0,n$. The points $z↓j$ vary smoothly, in
the sense that they lie on $m$ cubic polynomials $f↓k(j)$, where $z↓j=f↓k(j)$
when $x↓{k-1}≤j≤x↓k$. Here $0=x↓0<x↓1<\cdots<x↓m=n$, and the algorithm
minimizes $m$. It solves the problem in such a way that the $x↓k$ are integers,
but the piecewise cubics don't match up at the junction points.

The driver program incorporated in this version of the algorithm accepts data
as a sequence of $n+2$ characters that are either "+" or "0" or "-", specifying
differences or derivatives. For example, {\tt++00+00++00--0} means $z↑\prime↓0=+1$,
$a↓0=0$, $a↓1=1$, $a↓2=1$, $a↓3=1$, $a↓4=2$, $a↓5=2$, $a↓6=2$, $a↓7=3$, $a↓8=4$,
$a↓9=4$, $a↓{10}=4$, $a↓{11}=3$, $a↓{12}=2$, $z↑\prime↓{12}=0$, $n=12$.
;
require "⊂⊃⊂⊃" delimiters; "used for macros"
define # = ⊂;comment⊃; "used henceforth instead of quoted comments like this"
define nextline = ⊂('15&'12)⊃ # carriage-return and line-feed in print commands;
define thru = ⊂step 1 until⊃ # abbreviation for for clauses;
define DEBUGONLY = ⊂⊃ # changed to ⊂comment⊃ when not debugging;
define saf = ⊂safe⊃ # used when an array is believed to require no bounds checks;
DEBUGONLY redefine saf = ⊂⊃ # when debugging, belief turns to disbelief;
DEBUGONLY external procedure bail # the SAIL debugger in case of need;

comment The following variables are used for input and output by the algorithm;
integer n # size of input data;
define cap=100 # capacity of the arrays, must be $>n$;
real saf array aa[0:cap] # given data $a↓j$;
real z0p, znp # given derivatives at the endpoints;
real saf array z[0:cap] # the answers $z↓j$;

comment The following variables are used as indices during the calculations;
integer i,j,k,m # miscellaneous index variables;
comment Simplex procedure;

comment The following variables are used by the simplex procedure;
integer saf array nonbasic[1:4] # indices of the four current nonbasic variables;
real saf array t[0:4,-cap:cap] # if $v↓j$ is basic, we have
	$v↓j+\sum↓1↑4 t[i,j]v↓{nonbasic[i]}=t[0,j]$;
real saf array a,b,c,d[0:4] # $A+\sum↓1↑4 a[i]v↓{nonbasic[i]}=a[0]$, etc.;
real saf array bico[0:3] # $bico[j]={k-1\choose j}$;
real saf array f[0:cap] # $f↓k ≤ A{k-1\choose0}+B{k-1\choose1}+C{k-1\choose2}
	+D{k-1\choose3} ≤ f↓k+tol$;
boolean saf array basic[-cap:cap] # true if $v↓j$ is basic;
real tol # given tolerance;
integer kmax # maximum value of $k$;

procedure simplex;
begin comment This is a general procedure that determines $A$, $B$, $C$, $D$
satisfying the conditions stated above, for as large a value of $k$ as possible,
given \.{tol} and the $f$ values.
;

comment The following statements initialize the arrays;
a[0]←f[1]+tol; b[0]←f[2]-f[1]; c[0]←f[3]-f[2]; d[0]←f[4]-f[3]-c[0];
c[0]←c[0]-b[0]; d[0]←d[0]-c[0];
for j ← 1 thru 4 do
	begin nonbasic[j]←j; t[0,-j]←tol; t[1,-j]←t[2,-j]←t[3,-j]←t[4,-j]←0;
	t[j,-j]←1; t[j,0]←0; basic[-j]←true;
	t[1,j]←t[2,j]←t[3,j]←t[4,j]←0; basic[j]←false;
	end;
a[1]←1; a[2]←a[3]←a[4]←0;
b[1]←-1; b[2]←1; b[3]←b[4]←0;
c[1]←1; c[2]←-2; c[3]←1; c[4]←0;
d[1]←-1; d[2]←3; d[3]←-3; d[4]←1;
bico[0]←1; bico[1]←3; bico[2]←3; bico[3]←1; k←4;
while true do
	begin comment This loop attempts to increase $k$;
	comment First we add two new equations;
	k←k+1; if k>kmax then return;
	bico[3]←bico[3]+bico[2]; bico[2]←bico[2]+bico[1]; bico[1]←bico[1]+1;
	t[0,-k]←a[0]+bico[1]*b[0]+bico[2]*c[0]+bico[3]*d[0]-f[k];
	t[0,k]←tol-t[0,-k];
	for j←1 thru 4 do
		begin t[j,-k]←a[j]+bico[1]*b[j]+bico[2]*c[j]+bico[3]*d[j];
		t[j,k]←-t[j,-k];
		end;
	basic[k]←basic[-k]←true;
	comment Now we check if the new equations cause revisions to the tableau;
	if t[0,k]≤0 then j←k else if t[0,-k]≤0 then j←-k else continue;
	while true do
		begin comment Now $v↓j$ is the only source of infeasibility;
		real tmin # minimum coefficient;
		real rmin # minimum ratio while searching for pivot element;
		real r # test ratio;
		real x,y;
		integer jj,jjj,jp,ii # more indices;
		i←1; tmin←t[1,j];
		if t[2,j]<tmin then begin i←2; tmin←t[2,j] end;
		if t[3,j]<tmin then begin i←3; tmin←t[3,j] end;
		if t[4,j]<tmin then begin i←4; tmin←t[4,j] end;
		if tmin≥0 then
			begin if t[0,j]=0 then done else return;
			end;
		rmin←1000000000 # infinity (approximately);
		for jj←-k thru k do if t[i,jj]>0 then
			begin if (r←t[0,jj]/t[i,jj])≤rmin then
				begin if r=rmin then
					begin comment Avoid circling, by making sure
						that the rows of the full tableau
						are lexicographically positive;
					integer index # smallest column where rows
						jj and jjj differ;
					index←jjj; jjj↔jj;
					for ii←1 thru 4 do if ii≠i and
						nonbasic[ii]<index then
						begin if t[ii,jjj]*t[i,jj]≠
						t[ii,jj]*t[i,jjj] then
							index←nonbasic[ii];
						if t[ii,jjj]*t[i,jj]>
						t[ii,jj]*t[i,jjj] then jjj↔jj;
						end;
					if jjj>jj then jj←jjj # restore jj;
					end
				else	begin rmin←r;
					jjj←jj # currently smallest row;
					end;
				end;
			end;
		jj←nonbasic[i];
		comment Now pivot to make $v↓{jj}$ basic instead of $v↓{jjj}$;
		basic[jjj]←false; nonbasic[i]←jjj;
		x←1/t[i,jjj] # reciprocal of pivot element;
		for ii←0 thru 4 do
			begin t[ii,jj]←x*t[ii,jjj]; t[ii,jjj]←0;
			end;
		t[i,jj]←x;
		for jp←-k thru k do if (y←t[i,jp]) and basic[jp] then
			begin for ii←0 thru 4 do t[ii,jp]←t[ii,jp]-y*t[ii,jj];
			t[i,jp]←-y*x;
			end;
		for ii←0 thru 4 do if ii≠i then
			begin a[ii]←a[ii]-a[i]*t[ii,jj];
			b[ii]←b[ii]-b[i]*t[ii,jj];
			c[ii]←c[ii]-c[i]*t[ii,jj];
			d[ii]←d[ii]-d[i]*t[ii,jj];
			end;
		a[i]←-a[i]*x; b[i]←-b[i]*x; c[i]←-c[i]*x; d[i]←-d[i]*x;
		basic[jj]←true;
		if t[0,j]>0 then done;
		end;
	end;
end;
comment The main algorithm;

integer lastk # points up to here have been handled;
procedure main_algorithm;
begin comment The idea of the main algorithm is simply to call the simplex
procedure repeatedly, going as far with cubics as possible;
;
m←1; lastk←0; tol←1;
while true do
	begin comment Main interation, either returns or increases $m$;
	if n-lastk≤ 3 then
		begin print(nextline,"Final cubic is trivial."); done;
		end;
	for k←lastk thru n do f[k-lastk+1]←aa[k];
	kmax←n-lastk+1;
	simplex;
	k←k-1;
	print(nextline,"Cubic #",m," handles ",k," points.",nextline);
	z[lastk]←a[0]-.5;
	z[lastk+1]←a[0]+b[0]-.5;
	z[lastk+2]←a[0]+2*b[0]+c[0]-.5;
	j←lastk;
		begin comment Evaluation of the $z↓j$ with current cubic;
		real f,df,ddf,dddf # differences;
		f←z[j+2] # $f(j+2)$;
		df←b[0]+c[0] # $\Delta f(j+1)$;
		ddf←c[0] # $\Delta↑2 f(j)$;
		dddf←d[0] # $\Delta↑3 f$;
		j←j+3;
		while j<lastk+k do
			begin ddf←ddf+dddf # $\Delta↑2 f(j-2)$;
			df←df+ddf # $\Delta f(j-1)$;
			f←f+df # $f(j)$;
			z[j]←f;
			j←j+1;
			end;
		end;
	for j←lastk thru lastk+k-1 do print(z[j]);
	lastk←lastk+k;
	if lastk>n then done;
	m←m+1;
	end;
end;
comment The driver program;

string s # given string;
integer char # current input character;

setprint("unrnd4.tmp","B");
while true do
	begin print(nextline,"Input:"); s←inchwl;
	DEBUGONLY if s="b" then bail;
	if s=0 then done;
	setprint(null,"F"); print(s); setprint(null,"B") # echo the input;
	case (char←lop(s)) of begin
	["+"] z0p←+1;
	["0"] z0p←0;
	["-"] z0p←-1;
	else begin print("?"); z0p←0 end
	  end;
	n←0; aa[0]←0;
	while true do
		begin case (char←lop(s)) of begin
		["+"] aa[n+1]←aa[n]+1;
		["0"] aa[n+1]←aa[n];
		["-"] aa[n+1]←aa[n]-1;
		else begin print("?"); aa[n+1]←aa[n] end
		  end;
		if s=0 then done;
		n←n+1;
		if n≥cap then
			begin print(nextline,"Table size exceeded!"); done;
			end;
		end;
	znp←aa[n+1]-aa[n]; z[n]←aa[n];
	main_algorithm;
	end;
end